#Copyright (c) 2012 LUQMAN HAKIM BIN ABDUL HADI (csilhah@nus.edu.sg)
#
#Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files 
#(the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, 
#merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is 
#furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES 
#OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE 
#LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR 
#IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

#to run: perl BSA-genome-builder.pl <genome_build>
#example_1: perl BSA-genome-builder.pl hg19
#example_2: perl BSA-genome-builder.pl mm9

#!/usr/bin/perl
use strict;

my $genus = $ARGV[0];
my %accepted_chr = (
'chr1' => 1,
'chr2' => 1,
'chr3' => 1,
'chr4' => 1,
'chr5' => 1,
'chr6' => 1,
'chr7' => 1,
'chr8' => 1,
'chr9' => 1,
'chr10' => 1,
'chr11' => 1,
'chr12' => 1,
'chr13' => 1,
'chr14' => 1,
'chr15' => 1,
'chr16' => 1,
'chr17' => 1,
'chr18' => 1,
'chr19' => 1,
'chr20' => 1,
'chr21' => 1,
'chr22' => 1,
'chrX' => 1,
'chrY' => 1,
'chrM' => 1 );

print "getting the required database...";
system("wget ftp://hgdownload.cse.ucsc.edu/goldenPath/$genus/database/refFlat.txt.gz");
system("wget ftp://hgdownload.cse.ucsc.edu/goldenPath/$genus/database/nestedRepeats.txt.gz");
system("wget ftp://hgdownload.cse.ucsc.edu/goldenPath/$genus/database/cpgIslandExt.txt.gz");

print "unpacking...";
system("gunzip refFlat.txt.gz");
system("gunzip nestedRepeats.txt.gz");
system("gunzip cpgIslandExt.txt.gz");
print "done!\n";

system("mv refFlat.txt $genus.refFlat.txt");
system("mv nestedRepeats.txt $genus.nestedRepeats.txt");
system("mv cpgIslandExt.txt $genus.cpgIslandExt.txt");

print "Modifying the original database...\n";
system("perl -pi -e 's/,\t/\t/g' $genus.refFlat.txt"); #take out extra ','
system("sort -k 1,1 -k 2,2 $genus.refFlat.txt > $genus.refFlat.txt.sort");

my $gene_info_cpg_name = {};
my $gene_info_cpg_start = {};
my $gene_info_cpg_end = {};
my $gene_info_repeat_name = {};
my $gene_info_repeat_start = {};
my $gene_info_repeat_end = {};
my $gene_info_repeat_strand = {};
my $gene_info_repeat_class = {};
my $gene_info_repeat_family = {};

print "Building...";
open SRT, "$genus.refFlat.txt.sort";
open WRT, ">$genus.refFlat.txt.sort.mod";

my $prev_refseq = "";
my $count = 1;
while(<SRT>) {
	my ($tss, $prom, $prom3k, $first_intron, $endGene, @ex_strt, @ex_end);
	chomp $_;
	my ($gene_name, $refseq_id, $chr, $strand, $start, $end, 
		$start_cds, $end_cds, $number_of_exons, $exon_start, $exon_end) = split "\t", $_;

	if ($accepted_chr{$chr}!=1) {next;}

	if ($prev_refseq eq $refseq_id) {
		$prev_refseq = $refseq_id;
		$refseq_id = $refseq_id.".".$count; #add count for duplicates
		@ex_strt = split ",", $exon_start;
		@ex_end = split ",", $exon_end;
		if ($strand eq "+") { 
			$tss = $start;
			$prom = $tss - 2000;
			$prom3k = $tss - 3000;
			$first_intron = $ex_strt[1];
			$endGene = $end;
		}
		elsif ($strand eq "-") { 
			$tss = $end;
			$prom = $tss + 2000;
			$prom3k = $tss + 3000;
			$first_intron = $ex_end[$number_of_exons-2];
			$endGene = $start;
		}
	}

	else {
		$prev_refseq = $refseq_id;
		$count = 1;
		@ex_strt = split ",", $exon_start;
		@ex_end = split ",", $exon_end;
		if ($strand eq "+") { 
			$tss = $start;
			$prom = $tss - 2000;
			$prom3k = $tss - 3000;
			$first_intron = $ex_strt[1];
			$endGene = $end;
		}
		elsif ($strand eq "-") { 
			$tss = $end;
			$prom = $tss + 2000;
			$prom3k = $tss + 3000;
			$first_intron = $ex_end[$number_of_exons-2];
			$endGene = $start;
		}

	}

	my $w = "$gene_name\t$refseq_id\t$chr\t$strand\t$start\t$end\t$start_cds\t$end_cds\t$number_of_exons\t$exon_start\t$exon_end\t$tss\t$prom\t$prom3k\t$first_intron\t$endGene";
	print WRT $w."\n";
	$gene_info_cpg_name->{$chr}{$tss} = "";
	$gene_info_cpg_start->{$chr}{$tss} = "";
	$gene_info_cpg_end->{$chr}{$tss} = "";
	$gene_info_repeat_name->{$chr}{$tss} = "";
	$gene_info_repeat_start->{$chr}{$tss} = "";
	$gene_info_repeat_end->{$chr}{$tss} = "";
	$gene_info_repeat_strand->{$chr}{$tss} = "";
	$gene_info_repeat_class->{$chr}{$tss} = "";
	$gene_info_repeat_family->{$chr}{$tss} = "";
}
close SRT; close WRT;

#add CpG information
open MOD, "$genus.cpgIslandExt.txt";
while (my $cpgline = <MOD>) {
	my $tmp_line;
	my $tmp_dist = 10000000000000;	
	chomp $cpgline;
	
	my @arr = split "\t", $cpgline;

	if ($accepted_chr{$arr[1]}!=1) {next;}

	my $hsh = $gene_info_cpg_name->{$arr[1]};
	foreach my $tss (keys %$hsh) {
			#print $tss."\n";
			my $a = abs($tss - $arr[2]);
			if($a < $tmp_dist) { $tmp_line = $arr[1]."---".$tss."///"; $tmp_dist = $a}
			elsif($a == $tmp_dist) { $tmp_line .= $arr[1]."---".$tss."///";}
	}
	my @e_arr = split "///", $tmp_line;
	foreach(@e_arr) {	
		my @f_arr = split "---", $_;
		$gene_info_cpg_name->{$f_arr[0]}{$f_arr[1]} .= $arr[4].","; 
		$gene_info_cpg_start->{$f_arr[0]}{$f_arr[1]} .= $arr[2].","; 
		$gene_info_cpg_end->{$f_arr[0]}{$f_arr[1]} .= $arr[3].","; 
	}
}
close MOD;

#add Repeats information
open MOD, "$genus.nestedRepeats.txt";
while (my $repline = <MOD>) {
	my $tmp_line;
	my $tmp_dist = 10000000000000;
	chomp $repline;

	my @arr = split "\t", $repline;
	if ($accepted_chr{$arr[1]}!=1) {next;}

	my $hsh = $gene_info_repeat_name->{$arr[1]};
	foreach my $tss (keys %$hsh) {
			my $a = abs($tss - $arr[2]);
			if($a < $tmp_dist) { $tmp_line = $arr[1]."---".$tss."///"; $tmp_dist = $a}
			elsif($a == $tmp_dist) { $tmp_line .= $arr[1]."---".$tss."///";}
	}
	my @e_arr = split "///", $tmp_line;
	foreach(@e_arr) {	
		my @f_arr = split "---", $_;
		$gene_info_repeat_name->{$f_arr[0]}{$f_arr[1]} .= $arr[4].",";
		$gene_info_repeat_start->{$f_arr[0]}{$f_arr[1]} .= $arr[2].",";
		$gene_info_repeat_end->{$f_arr[0]}{$f_arr[1]} .= $arr[3].",";
		$gene_info_repeat_strand->{$f_arr[0]}{$f_arr[1]} .= $arr[6].",";
		$gene_info_repeat_class->{$f_arr[0]}{$f_arr[1]} .= $arr[15].",";
		$gene_info_repeat_family->{$f_arr[0]}{$f_arr[1]} .= $arr[16].",";
	}
}
close MOD;

open GENE, "$genus.refFlat.txt.sort.mod";
open WRT, ">" , "$genus.refFlat.w_cpg.txt";
while(my $line = <GENE>) {
	chomp $line;
	my @arr = split "\t", $line;
	print WRT $line."\t".$gene_info_cpg_name->{$arr[2]}{$arr[11]}."\t".$gene_info_cpg_start->{$arr[2]}{$arr[11]}."\t".$gene_info_cpg_end->{$arr[2]}{$arr[11]}."\t".$gene_info_repeat_name->{$arr[2]}{$arr[11]}."\t".$gene_info_repeat_start->{$arr[2]}{$arr[11]}."\t".$gene_info_repeat_end->{$arr[2]}{$arr[11]}."\t".$gene_info_repeat_strand->{$arr[2]}{$arr[11]}."\t".$gene_info_repeat_class->{$arr[2]}{$arr[11]}."\t".$gene_info_repeat_family->{$arr[2]}{$arr[11]}."\n";
}
open HDR, ">header.txt";  
print HDR "Gene_name\tGene_ID\tchr\tstrand\tstart\tend\tStart-cds\tEnd-cds\tnb_exon\texon_start\texon_end\tTSS\tprom\tprom3k\tfIntron\tendGene\tCpG_name\tCpG_start\tCpG_end\tRepeat_name\tRepeat_start\tRepeat_end\tRepeat_strand\tRepeat_class\tRepeat_family\n";
close HDR;
system("perl -pi -e 's/,\t/\t/g' $genus.refFlat.w_cpg.txt");
system("sort -k 3,3 -k 12,12n $genus.refFlat.w_cpg.txt > $genus.refFlat.w_cpg.txt.tmp");
system("mv $genus.refFlat.w_cpg.txt.tmp $genus.refFlat.w_cpg.txt");

system("cat header.txt $genus.refFlat.w_cpg.txt > $genus.refFlat.w_cpg.txt.tmp");
system("mv $genus.refFlat.w_cpg.txt.tmp $genus.refFlat.w_cpg.txt");

print "done!\n";